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PREFACE 


Many  problems  in  modern  control  engineering  involve  the  numeri¬ 
cal  solution  of  nonlinear  differential  equations  subject  to  two- 
point  boundary  conditions.  Hie  purpose  of  this  Memorandum  is  to 
show  how  nonlinear  extrapolation  techniques  can  be  used  to  convert 
a  first-order  successive  approximation  scheme  into  a  second-order 
scheme  for  the  solution  of  such  problems.  Results  of  a  numerical 
experiment  are  presented. 
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SUMMARY 


It  is  suggested  that  the  convergence  properties  of  the  usual 
Picard  successive  approximation  scheme  may  be  improved  through  use  of 
nonlinear  extrapolation  techniques.  A  numerical  example  is  provided. 
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I.  INTRODUCTION 


An  interesting  class  of  nonlinear  two-point  boundary  value 
problems  is  described  by  the  following  equations 


u"  =  f(u),  u(0)  =  u(l)  =  0. 


(1) 


In  the  scalar  case,  the  equation  can  be  solved  by  quadrature;  in 
the  vector  case,  we  must  use  numerical  techniques.  Let  us  restrain 
our  attention  to  the  scalar  case  in  order  to  present  a  new  approach 
unhindered  by  analytic  details.  One  way  to  attack  these  problems 
is  to  use  Newton's  method  and  quasilinearization. ^  Such  an  ap¬ 
proach  requires  knowledge  of  the  partial  derivative  cjf/du  for  the 
determination  of  the  (n+l)st  approximation  to  the  solution,  given 
an  n  approximation,  u^(x) .  The  recurrence  relations  are  the  linear 
equations 


u"  .  =  f(u  )  +  (u  -  u  )  s  - 
n+1  v  n'  n+1  n'  ^ 


(2) 


un+l(°>  =  W1*  =  °* 


They  can  be  solved  by  producing  a  particular  solution  and  two  in¬ 
dependent  solutions  of  the  homogeneous  equations  numerically  and 
determining  the  constant  multipliers  of  the  homogeneous  solutions 
so  as  to  satisfy  the  boundary  conditions.  The  sequence  of  functions 
u  (x),  u^x),  u^x),  ...  converges  quadratically  to  the  solution 
u(x),  0  =  x  =  1,  under  appropriate  conditions  on  the  function  f(u); 


cf.  Ref.  1. 
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s  t 

Picard's  method  for  finding  the  (n+1)  approximation  requires 
the  solution  of  the  linear  system 


u"  .  =  f(u  ), 
n+1  v  n  " 


un+l(0)  =  un+l(l)  =  °* 


(3) 


s  t 

It  is  a  first  order  method,  so  that  if  the  (n+l)°  approximation  is 


u  =  u  +  5  , 
n+1  n’ 


(M 


vhere  8^  is  the  discrepancy,  then 


6n+lU~K!  5n  II* 


(5) 


If  K  <  1,  the  sequence  of  approximating  functions  converges  linearly 
or  "geometrically." 


II.  SKETCH  OF  METHOD 

Let  us  assume  that  we  have  a  current  approximation  UQ(x), 

0  =  x  =  1.  We  would  like  to  find  a  new  approximation  using  a  varia¬ 
tion  of  Picard's  method.  The  method  which  we  propose  has  the  ad¬ 
vantage  of  quadratic  convergence  without  requiring  any  partial 
differentiation  which  could  require  a  major  calculation  in  a  larger 
problem.  Let  us  compute  the  functions  u^(x)  and  Ug(x)  by  means  of 
the  equations 

u"  =  f(UQ),  u1(0)  =  u^l)  =  0, 

u2  =  f(ul)>  ^(O)  =  u2(l)  =  0. 


(6) 
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An  initial  approximation  was  computed  with  u(0)  =  0,  u'(0)  =  -  1, 
and  the  differential  Ea.  (1).  The  integration  step  length  was  .01. 
The  results  of  this  experiment  are  given  in  the  table,  together 
with  the  correct  solution  u(x). 


X 

UQ(x) 

Ux(x) 

U4(x) 

u(x) 

0.0 

0.000000 

0.000000 

0.000000 

0.000000 

0.1 

-0.095159 

-0.041328 

-0.04l487 

-0.o4i436 

0.2 

-0.181212 

-0.073067 

-0.073337 

-0.073269 

0.3 

-O.258912 

-0.095566 

-0.095861 

-0.095800 

0.4 

-0.328885 

-0.10906l 

-0.109274 

-0.109238 

0.5 

-0.391654 

-O.H3677 

-0.113708 

-0.113704 

0.6 

-0.447657 

-0.109433 

-0.109209 

-0.109238 

0.7 

-0.497265 

-0.096231 

-0.095746 

-0.095800 

0.8 

-0.540786 

-0.073855 

-0.073205 

-0.073269 

0.9 

-0.578481 

-0.041950 

-0.041387 

-o.o4i436 

1.0 

-0.610565 

0.000000 

0.000000 

0.000000 

Note  that  a  single  iteration  gave  a  very  great  improvement  over  the 
initial  approximation.  The  time  required  to  do  these  four  itera¬ 
tions  on  the  IBM  7090  computer  was  about  20  seconds. 

The  FORTRAN  II  programs  are  listed  below. 


REFERENCES 


1.  Kalaba,  R.  E.,  "On  Nonlinear  Differential  Equations,  the 

Maximum  Operation  and  Monotone  Convergence,"  J.  of  Math  and 
Mech.,  Vol.  8  (1959),  pp.  519-574. 

2.  Shanks,  Daniel,  "Non-linear  Transformations  of  Divergent  and 

Slowly  Convergent  Sequences,"  J.  Math,  and  Phys.,  Vol.  34 
(1955),  PP.  1-42. 


U  J  U  KJ  U  KJ  U  KJ  J 


5 


*  LIST 

CPICARD  PICARD'S  METHOD  ♦  2-PT .  B.V. PROBLEMS 

C 

COMMON  T  *NMAX  .KMAX  »N1 .NPRNT  .UZERO.U. I  FLAG  » DELTA »UP REV 
DIMENSION  T( luO) .UZEPOI2) .UI2G1) 

DIMENSION  PI  I  201 } .P2 I  201 ) 

DIMENSION  U1I201) 

C 

1  READluOO.NMAX.KMAX.Nl .NPRNT 
PRINT910 

PRINT90G.NMAX.KMAX.N1 .NPRNT 
RFADl^OI .DELTA. (U2ER0I 1) .1*1.2 » 

PR1NT901. DELTA. (UZEROi 1 ) . 1=1.2) 

C 

K*  0 

PRINT902.K 
CALL  START 

K  ITERATIONS,  EACH  a  DOUBLE  APPLICATION  OF  PICARD'S  METHOD 

DO  100  K=  1  * KMAX 
PRINT902.K 
I  FLAG* 1 

APPLICATION  1 
T ( 2 ) “0  « 

T ( 3 ) *DELT  A 
TC4)»0- 
T ( 5 ) *0 • 

UPREV«U( 1 ) 

(_au_L  INT(T.  2.  Nl,  0.0.0. 0.0.0) 

PI < 1  )  *T I  4 ) 

DO  2  N-2.NMAX 
UPKt V*U ( N— 1 ) 

CALL  I NTM 

2  P 1  (  N  )  =  T  (  4  ) 

Bl—Pl  (NMAX) 

APPLICATION  2 
T ( 2 ) *C  # 
r<4)«o. 

TC5»«0. 

U1(1)=P1(1) 

UPREV=U1 ( 1 ’ 

CALL  INTI T. 2, Nl ,0.0. 0.0. 0.0  ) 

P2 1 1 ) *T I  4 ) 

C 

DO  3  N=2»NMAX 
UPRtVUl  (N-i  ) 

CAI  L  I  NTM 

U1(N)*P1(N)  +  B1*T I  2 ) 

3  P2(N)«TI4) 


n  n  n 
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C 

52  E_P  2 ( NM  A  X ) 

NEW  APPROXIMATION 

PR INT903 
U( 1 ) =w. 

T ( 2  J  =o  . 

PR INT904, T ( 2)  .Ut 1 ) 

N*  1 
C 

5  DO  53  M*1.NPRNT 
N  =  N+1 

T(2)*T(2)+DELTA 
U  2 :t  P  2  (  N )  +  B2*T(2) 

I F ( N-NMAX )51»52*52 

51  U ( N ) *  U  2  -  (U1(N)-U2>**2  /  < U ( N ) -2 .*U1 I N ) +U2  ) 

GO  TO  53 

52  U ( NMA X  )  *0  # 

53  CONTINUE 
C 

PR  I NT9C4»  T ( 2 )  »U(N) 

I F ( N-NMAX ) 5*1^0*100 
C 

100  CONTINUE 
C 

GO  TO  1 
C 

1000  FORMAT (6112) 

1001  FORMAT  <  6E1 2  •  8  ) 

900  FORMAT (6120) 

901  FORMAT (6E20. 6  ) 

902  FORMAT  1 10H1  ITERAT ION»  13) 

903  FORMAT ( 1H018X1HT  *26X4HU( T )  /> 

904  FORMAT ( 10XF10. 4,  E30.8) 

910  FORMAT ( 3 7H 1 P 1 C ARD ' S  METHOD  ♦  2-PT .  B . V . PROBLEMS  ) 
END 
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*  LIST 

SUBROUTINE  START 

COMMON  T.NMAX  »KMAX.N1 »NPRNT  *UZERO .U . I F LAG » DELTA . UPREV 
DIMENSION  TC luO) »UZERO(2) ,U<201) 

ACTUAL  EQS. 

IFLAG-2 
T(2)-C. 

T ( 3 ) *DELT  A 
TI4)»llZFRO(  1  \ 

T ( 5 ) *UZERO  <  2 ) 

CALL  INT(T,2.N1.0.0»0.0.0»0) 

N-l 

U(1I«TI4) 

PRINT903 

PR INT904»  T ( 2  I »U(1 ) 

C 

5  DO  6  M«1,NPRNT 
CALL  I NTM 
N-N+l 

6  U  <  N  )  ■  T  ( 4  ) 

C 

PRINT904.T ( 2) »U(N) 

I F ( N-NMAX ) 5  »7  *7 

7  RETURN 
C 

903  FORMAT (1H018X1HT  »  26X4HU  t  T  >  /) 

904  FORMAT ( 10XF10»4*  E30.8) 

END 
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LIST 

SUBROUTINE  DAUX 

COMMON  T.NMAX .KMAX.Nl .NPRNT .UZERO . U  *  I  FLAG . DELT A .UPREV 
DIMENSION  T  < 1 U 0 ) *UZ ERO { 2 ) ♦ U ( 201 ) 

GO  TO  ( 1*2 ) » I  FLAG 

PTCARD'S  EOS. 

T  <  6 ) =T ( 5  ) 

T ( 7 ) =EXPF (UPREV) 

RETURN 


ACTUAL  EQS. 
T(6»«T(5) 

T ( 7 ) *EXPF ( T  t  4 ) ) 

RETURN 

END 
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DATA 

201 

.OOS 


0 
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